#analyze_id.pl
#author: Jeremy Cothran (jeremy.cothran@gmail.com)
#date: March 28, 2008

use strict;

my $summary_file = 'summary_id.txt';

if (@ARGV < 1) { print "usage like : perl analyze_id.pl query_out.txt\n"; exit 0; }

my ($blast_file) = @ARGV;


###############################################

my $cutoff_percent = 94;
my $mismatch_tolerance = 0;  #number of bases allowed to mismatch


open (FILE_IN,"$blast_file");

my %hash_seq = ();
my $r_hash_seq = \%hash_seq;

my ($junk0,$seq_count,$seq_interest,$length_seq_interest,$seq_identity);

my $line;
while ($line = <FILE_IN>) {

	if ($line =~ /^analysis/) {

		($junk0,$seq_count,$seq_interest) = split(/\s+/, $line);
		$length_seq_interest = length($seq_interest);
		
		next;
	}

	if ($line =~ /^identity/) {

		$seq_identity = substr($line,9,67);
		
		next;
	}


	if ($line =~ /^tmpseq_0/) {

		my ($junk0,$acc_id,$percent,$length_match) = split(/\s+/, $line);

		$hash_seq{'all'}{seq}{$seq_interest}{count} = $seq_count;
		$hash_seq{'all'}{seq}{$seq_interest}{identity} = $seq_identity;
		$hash_seq{'all'}{seq}{$seq_interest}{percent} = $percent;
		$hash_seq{'all'}{seq}{$seq_interest}{length_match} = $length_match;

		my $mismatch = $length_seq_interest - $length_match;

		if ($mismatch <= $mismatch_tolerance) {
			if ($percent >= $cutoff_percent){
				$hash_seq{'all'}{seq}{$seq_interest}{quality} = 'good';				
			}
		}
		else { $hash_seq{'all'}{seq}{$seq_interest}{quality} = 'bad'; }
			
		next;

	} #if $line

} #while $line

close (FILE_IN);

###############################################

#summary


open (FILE_OUT, ">$summary_file");
print FILE_OUT "#cutoff_percent=$cutoff_percent mismatch_tolerance=$mismatch_tolerance\n";
print FILE_OUT "#quality(+++=good,---=bad)\tcount\tidentity\tmatch_percent\tmatch_align\n";

#order by sample count
my @seq = reverse sort { $hash_seq{'all'}{seq}{$a}{count} <=> $hash_seq{'all'}{seq}{$b}{count} } keys %{$r_hash_seq->{'all'}{seq}};

foreach my $seq (@seq) {
	
	my $length_seq = length($seq);
	my $length_match = $hash_seq{'all'}{seq}{$seq}{length_match};
	
	if ($hash_seq{'all'}{seq}{$seq}{quality} eq 'good') {
		print FILE_OUT "+++\t$hash_seq{'all'}{seq}{$seq}{count}\t$hash_seq{'all'}{seq}{$seq}{identity}\t$hash_seq{'all'}{seq}{$seq}{percent}\n";
	}
	else {
		print FILE_OUT "---\t$hash_seq{'all'}{seq}{$seq}{count}\t$hash_seq{'all'}{seq}{$seq}{identity}\t$hash_seq{'all'}{seq}{$seq}{percent}\t$length_match/$length_seq\n";
	}
		
}

close (FILE_OUT);

##############################################################################################
###############################################

my @unique_files = glob 'found/unique_found*';

while (my $sub_file = shift @unique_files) {

my $sub_report = 'report_'.substr($sub_file,6);

$hash_seq{'report_good_97_100'}{seq} = ();
$hash_seq{'report_good_94_97'}{seq} = ();
$hash_seq{'report_bad'}{seq} = ();

my $total_sample_count = 0;


###############################################

#hash report populate

open (FILE_IN,"$sub_file");

my $line;
while ($line = <FILE_IN>) {

	my ($seq_count,$seq_interest) = split(/\s+/, $line);

	my $seq_identity = $hash_seq{'all'}{seq}{$seq_interest}{identity};
	if (!($seq_identity)) { $seq_identity = 'Unprocessed'; } #sequences still in process and not part of final summary yet

	my $quality = $hash_seq{'all'}{seq}{$seq_interest}{quality};
	my $percent = $hash_seq{'all'}{seq}{$seq_interest}{percent};
	my $length_match = $hash_seq{'all'}{seq}{$seq_interest}{length_match};

	if ($quality eq 'good') {
		if ($percent >= 97) {
			$hash_seq{'report_good_97_100'}{seq}{$seq_identity}{count} += $seq_count;
		}
		elsif ($percent >= 94) {
			$hash_seq{'report_good_94_97'}{seq}{$seq_identity}{count} += $seq_count;		
		}
	}
	else {

		$hash_seq{'report_bad'}{seq}{$seq_identity}{count} += $seq_count;
	
	}
	
	$total_sample_count += $seq_count;
	
} #while $line

close (FILE_IN);

###############################################

#hash sort/report

open (FILE_OUT, ">found/$sub_report");
print FILE_OUT "#$sub_file\n#sample_count\tsample_percent(of $total_sample_count samples)\tidentity\n";

my @array_report_band = ('report_good_97_100','report_good_94_97','report_bad');

foreach my $report_band (@array_report_band) {

print FILE_OUT "#$report_band\n";

#order by sample count
my @seq = reverse sort { $hash_seq{$report_band}{seq}{$a}{count} <=> $hash_seq{$report_band}{seq}{$b}{count} } keys %{$r_hash_seq->{$report_band}{seq}};
foreach my $seq (@seq) {

	my $sample_count = $hash_seq{$report_band}{seq}{$seq}{count};
	
	my $sample_percent = sprintf("%.3f",$sample_count/$total_sample_count*100);
	
	print FILE_OUT "$sample_count\t$sample_percent\t$seq\n";
		
}

} #foreach $report_band

close (FILE_OUT);

} #while $unique_files

exit 0;
